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We calculate the energies of quarkonium bound states in the presence of a medium of nonzero 
isospin density using lattice QCD. The medium, created using a canonical (fixed isospin charge) 
approach, induces a reduction of the quarkonium energies. As the isospin density increases, the 
energy shifts first increase linearly and then saturate at an isospin density where a qualitative 
change in the behaviour of the energy density of the medium has been previously observed. This 
provides further support to the conjecture that the medium undergoes a transition from a pion 
gas to a Bose-Einstein condensed phase at this point. The reduction of the quarkonium energies 
becomes more pronounced as the heavy-quark mass is decreased, similar to the behaviour seen in 
two-colour QCD at non-zero quark chemical potential. In the process of our analysis, the T-tt and 
?7(,-7r scattering phase shifts are determined at low momentum and found to be small. 



I. INTRODUCTION 

An important probe of exotic phases of QCD matter 
is the way in which heavy quarkonium propagation is 
modified by by the presence of that matter. The heavy 
quarks can in some sense be viewed as separable from the 
medium which is predominantly composed of light quark 
and gluonic degrees of freedom. At non-zero tempera- 
ture, the suppression of the propagation of J/ij} parti- 
cles is a key signature for the formation of a quark-gluon 
plasma [1]. This suppression has been observed for char- 
monium in various experiments at SPS, RHIC and the 
LHC and recently in the T spectrum at the LHC [5]. 
Quarkonium propagation is naturally also expected to 
be a sensitive probe of other changes of phase such as 
those that occur at high density or large isospin density. 

Since the effects of QCD matter on quarkonia are es- 
sentially non-perturbative in origin, a systematic eval- 
uation requires input from lattice QCD. At some level, 
these effects can be distilled to a change in the potential 
between the quark-anti-quark {QQ) pair that binds them 
into quarkonium. At non-zero temperature but zero den- 
sity, this has been studied extensively using lattice QCD 
(see Ref. [5] for a recent overview) where strong screening 
effects are seen near the deconfinement scale. Significant 
effects are also seen in investigations of the properties 
of charmonium and bottomonium spectral functions at 
non-zero temperature (see OH]). 

Modifications of the potential or quarkonium proper- 
ties will also occur for non-zero density. Ref. [5_ has in- 
vestigated the static potential in the presence of a gas 
of pions and below we briefly address how the medium 
affects the binding of quarkonium through solving the 
Schrodingier equation for the modified potential. As 
the main focus of this work, however, we explore the 
effects of isospin charge density on quarkonium bound 
state energies more directly by using lattice NRQCD 
(non-relativistic QCD) to compute quarkonium correla- 
tion functions in the presence of a medium of varying 
isospin chemical potential. At low isospin densities, and 



correspondingly low chemical potentials, we find that the 
ground state energy of the quarkonium systems decreases 
with increasing density, showing qualitative agreement 
between the potential model calculation and the QCD 
calculation. However, at an effective isospin chemical 
potential /// ^ /i/,pcak = 1.3 ttt-tt (where previous calcu- 
lations of the energy density of the isospin medium have 
suggested a transition to a Bose-Einstein condensed state 
[5] in line with theoretical expectations [7]), the effect of 
the medium on the quarkonium energy saturates to a 
constant shift. At still larger chemical potentials, the 
determination of the energy shift becomes statistically 
noisy. The marked change in behaviour provides addi- 
tional support to the notion of the change in phase of 
the system. 

Our study is presented as follows. In Section |lTj we 
construct an isospin density dependent potential model 
of quarkonium states, and calculate the expected shift 
of the quarkonium energies in the medium. Section |III| 
presents the methodology of the lattice QCD calculation 
of these energy shifts, with results presented in Sections 
IIVI and |Vl Wc conclude with a discussion of the results 
and possible extensions to the current work. 



II. MODIFICATION OF THE STATIC 
QUARK-ANTI-QUARK POTENTIAL 

Here we briefly review the results obtained in Ref. [S] 
on the screening of the static potential by a pion gas and 
ascertain the expected effects of these modiflcations on 
quarkonium spectroscopy. Ref. presented a sophisti- 
cated calculation of the static quark-anti-quark potential 
in a pionic medium as a function of the separation of the 
QQ pair, r, and as a function of the isospin charge, n, of 
the medium. This required measurement of Wilson loops 
of various extents in space and in time and to obtain sig- 
nals for large loops, different levels of HYP (hypercubic) 
smearing [8] were used in overlapping regions of r, com- 
plicating the analysis. As will be discussed below, con- 
structing appropriate ratios of correlation functions was 
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also critical in order to obtain statistically clean mea- 
surements. The central results of this work were that the 
potential is screened by the presence of the medium and 
that this screening effect is small. For the relatively low 
pion densities investigated in Ref. 5J , the dominant effect 
corresponded to a change in the potential in the linearly 
rising region that was approximately linearly dependent 
on both r and n. This form, dV{pi^r) = a pi r, cor- 
responds to the physical expectation of a gas of weakly- 
interacting pions permeating a flux-tube of constant ra- 
dius between the static quark and anti-quark, and is a 
picture in which the appearance of the isospin density. 
Pi, is natural. Performing a correlated fit to the results 
presented in Ref. [5 using this form, we are able to de- 
scribe the data well, as shown in Fig. [T] and find that 
a = -9.6(1.6) MeV fm^. 

To estimate the effects on quarkonium spectroscopy. 



we use the Cornell potential Vcorncii(f') — — (4/3)as/r + 
K r with as = 0.24 and ^/K = 468 MeV (values fixed 
in vacuum from Ref. ^9j) and augment it with the small 
screening shift discussed above. We then solve the radi- 
ally symmetric Schrodinger equation numerically for an- 
gular momentum i and reduced mass mred = m/2 (where 
m is the heavy-quark mass), 
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£{£ + 1) 



+Vcornou(?') + 5V{pi,r) ui{r) = E ui{r) , 



(1) 



to establish the wave functions and eigenstate energies 
for the various quantum numbers. The energy shift is 
then defined simply as the difference of the resulting en- 
ergy from that where 5V{pi,r) is omitted. We calculate 
this shift for both the 15' and IP states and various dif- 
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FIG. 2: Shifts in the energies of the 15* (upper) and IP (lower) 
states in quarkonium as a function of the isospin density. Re- 
sults are shown for four different values of the heavy-quark 
mass with the uncertainty shown only for the mass closest to 
the physical bottom-quark mass, m = 4.676 GeV (uncertain- 
ties for the other masses are of similar size). 



ferent values of the heavy-quark mass as shown in Fig. [2] 
We note that one could use only the Cornell potential 
to determine the wave functions and include the addi- 
tional small shift from the screening as a perturbation, 
calculating 

5E{pj) ^ Jdr uf^*{r)SV{pi,r)uf\r) , (2) 

where uf\r) are the solutions to Eq. when SV{pi,r) 
is omitted. As the shift is small, we expect this will give 
consistent results. 

Since P-wave states are more extended in size, they 
probe regions of the potential where the shift is larger 
and consequently we find that the energy shift is larger 
for these states than for the S-wave states. The effect 
also increases as the heavy-quark mass decreases, again 
because of the larger size of the lighter systems. 



In the following, we determine quarkonium eigen- 
energies in (NR)QCD at non-zero isospin density and 
investigate to what extent they are predicted by the po- 
tential model described above based on a screening pion 
gas. 

III. LATTICE METHODOLOGY 

A. Lattice details 

In this study, we make use of anisotropic gauge 
configurations generated by the Hadron Spectrum and 
NPLQCD collaborations. The full details of the action 
and algorithms used to generate the configurations are 
discussed in the original works, Refs. [TUl HI]; here we 
summarise the salient features of the configurations and 
the measurements that we perform. A tree-level, tadpole- 
improved gauge action |12| , and nf — 2 + 1 flavour clover 
fermion action |13j are used. Two levels of stout smear- 
ing [H] with weight p — 0.14 are applied in spatial di- 
rections only in order to preserve the ultra-locality of the 
action in the temporal direction. The gauge action is 
constructed without a 1 x 2 rectangle in the time direc- 
tion for the same reason. In this study, we make use of 
a single spatial lattice spacing, = 0.1227(8) fm [TT] 
and have a renormalised anisotropy of f = Ug/at = 3.5, 
where at is the temporal lattice spacing. We also work at 
a single value of the light-quark mass for this exploratory 
investigation and use a strange-quark mass that is close 
to its physical value; these values correspond to a pion 
mass of ~ 390 MeV and a kaon mass of rriK ~ 540 
MeV. For these parameters, we investigate three different 
ensembles, corresponding to different physical volumes 
and temporal extents as shown in Table |l] The differ- 
ent physical volumes allow us to access a large range of 
isospin densities in our study, and the different tempo- 
ral extents provide control of thermal effects as discussed 
in Ref. [B]. On these gauge configurations we calculate 
correlation functions involving light quarks and use the 
colourwave propagator basis introduced in Ref. [6l, fix- 
ing to Coulomb gauge and using plane- wave sources and 
sinks for a range of low momenta (A'mom in total on each 
ensemble, see Table |l]). For each case, we calculate light- 
quark propagators on A"cfg configurations from Agrc time- 
slices, equally spaced throughout the temporal extent. 
Details of the NRQCD heavy quark propagator calcula- 
tions are discussed below. 
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TABLE I: Details of the ensembles and measurements used 
in this work, uq^ is defined as the fourth root of the spatial 
plaquette. 
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B. Multi-pion lattice correlators 

In order to produce the medium that will modify the 
propagation of the quarkonium states, we use the canoni- 
cal approach of constructing many-pion correlation func- 
tions that is described in detail in Ref. [5], using meth- 
ods developed there and in earlier works [T5HT9] . As 
discussed therein, correlators of a fixed isospin charge, 
n = X^iLi^-j' total momentum, Pj, making use of 
N sources, are given by 

\i=l \xi,x' 



xH E^''''^''^"(y^'0)^5d(y,'0)| ),(3) 

where P f — P/ii and, for momentum conservation, 

J2iLiiPi ~ P2) = Sj^i Pfj ■ In what follows, we will set 
P f = but the p/. and 2 take various values subject 
to these constraints; different choices of the momenta 
defining the interpolating operators will have different 
overlap onto the eigenstates of the chosen P / but provide 
additional statistical resolution in the determination of 
the energy of the system. 

To construct these correlation functions, we work in 
Coulomb gauge and compute light-quark colourwave 
propagators 
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FIG. 3: riij and T correlators (upper) and effective energies 
(lower) on the 20^ x 256 ensemble, for asin — 2.6. 



where 

S,/<i(x, t; p', 0) = ^ e'P'yS,/<i(x, t; y, 0) 
y 

is a solution of the lattice Dirac equation: 

D{y, t- X, t)5„/d(x, t- p', 0) = e^P'^'^f^o . 

x,t 

The contractions implicit in Eq. ([s]) can be written in 
terms of a matrix A, the 12 x 12 sub-blocks of which are 
given by 

Ak, [t) ^Y.^ (Pi ' P) (-P2' P/. - P) ' (5) 



and can be computed efficiently using the methods of 



Ref. 0. The diflFerent C„ 



for a given n occur in 



complicated combinations in this expansion, however we 
are explicitly only interested in the energies of the system, 
so the particulars of the combination are irrelevant. 

These correlators have been studied in detail in pre- 
vious work [5] and we do not present them again here. 
As investigated in detail in Ref. [6], many-pion corre- 
lations contain thermal contributions in which parts of 
the system propagate around the temporal boundary. In 
our choice of fitting ranges in the analysis presented be- 
low, we are careful to remain away from the regions in 
Euclidean time that are contaminated by either excited 
states or by these thermal effects. 



where k, i label the source and sink, and the dependence 
on pj, P2, and P/ is suppressed. The correlators above 
can be extracted by noting that combinations of the 
C'ni, for a given n ~ X^iLi are the coefficients 



of the expansion of 

det[l -f Ai] = 1 + ACi, + \^C2^ + ■■■ + Ai2^Ci2A.. ,(6) 



C. NRQCD for quarkonium correlators 

To implement the heavy quarks in our quarkonia sys- 
tems, we use a lattice discretisation of non-rclativistic 
QCD (NRQCD). Since our light quark and gluon de- 
grees of freedom are defined on an anisotropic lattice. 
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FIG. 4: /lb correlator (upper) and effective energy (lower) on 
the 16^ X 128 ensemble, for aaVa = 2.6. 



we require lattice NRQCD [201 [21] formulated on an 
anisotropic lattice as first set out in Ref. [52] . As the non- 
relativistic nature of the theory already separates space 
and time, using a temporal lattice spacing that differs 
from the spatial lattice spacing is a very natural choice for 
NRQCD. Anisotropic lattice NRQCD has been used for 
example to calculate the spectrum of quarkonium hybrid 
states [531 [23] , and recently also to study quarkonium at 
non-zero temperature [3J [5S| . 

The Euclidean action for the heavy quark field, "0, can 
be written as 



al V V^(x, t) [V(x, t) ~ K{t) V'(x, t - at)] , (7) 



where K[t) is the operator that evolves the heavy-quark 
Green function forward one step in time. Here we use 
the form 



m = 1 - 



at 5H\ 



1 - 



2 J \ 2n 

0.tHo\t-at 



1 - 



2n 



1 - 



U^t-at) 
at SH\t-at 



, (8) 



where Uq are the temporal gauge links. In this expres- 
sion, 



Hn = 



A(2) 



(9) 



is the order-w^ term in the NRQCD velocity expansion, 
and SH is a correction term given by 



SH = -ci ^ ^ J + C2 ^ f V • E - E • V 



8m3 



8m2 



-C3 ^ cr ■ f V X E - E X V) - C4 cr • B 

8m^ \ / Zm 



-C5 



24m 



C6 



at 



(A(2)) 



16n 



(10) 



(the notation is as in Ref. [26]). The operators with coef- 
ficients Ci through C4 are the relativistic corrections of or- 
der v"^, the operator with coefficient C5 removes the order- 
discretization error of Hq, and the operator with cg 
removes the leading order-at error in the time evolution. 
We use the tree-level values of the matching coefficients, 
Ci = 1. For tadpole improvement [27] of the derivatives 
and field strengths, we set uqs equal to the 4th root of 
the spatial plaquette (see Table [l]), and set uot ~ 1- 

To avoid instabilities in the time evolution with the 
operator in Eq. ([8]) , the parameter n must be chosen such 
that max[aii/o/(2ri)] < 2 [21]. On an anisotropic lattice, 
this requires 



n > 3at/{2alm) = 3/(2^0^771) 



(11) 



(interactions with gluons weaken this requirement 
slightly [H]). In this work, we set the bare heavy-quark 
mass to a^TO = 2.6 (which is near the b quark mass) as 
well as to the lower values astn = 1.8,1.4,1.1. Because 
we have ^ = 3.5, a stability parameter of n = 1 is suffi- 
cient in all cases. 

When using NRQCD, all quarkonium energies are 
shifted by an unknown constant (which is approximately 
equal to two times the heavy-quark mass). This shift is 
state-independent and cancels in energy splittings as well 
as in differences between energies extracted at zero and 
non-zero isospin density. For the purpose of tuning the 
heavy-quark mass, we measure the kinetic masses of the 
rji, and T states, defined as 



at Mkin = 



iasP)ye - [atEjp) - atEjof 
2 [atE{p) - atEiO)] 



(12) 



with one unit of lattice momentum, |p| — 2tt/L. The 
spin-averaged values of the IS* kinetic masses, Afkin = 
(3Af]J„ + M^^J/A computed on the 16^ x 128 ensemble 
(at p/ = 0), are given in Table [iTj 

In the main calculations of this work, we use zero- 
momentum smeared quarkonium interpolating fields of 
the form 

ouit) = E E (y'' ^)r(y - y')^iy^ t) (is) 

y' y 
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0.8732(49) 
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1.0138(20) 
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1.8 


1.2225(15) 


6.881(46) 


2.6 


1.6590(12) 


9.337(62) 



TABLE II: Spin-averaged quarkonium kinetic masses on the 
16^ X 128 ensemble. 



at the sink and 

0-66(0) = E^^(O'0)rWV'(x,0) (14) 

X 

at the source. Here, x is the heavy anti-quark field 
and r(r) is the smearing function, which is a 2 x 2 
matrix in spinor space. Note that antiquark propaga- 
tors can be obtained from quark propagators through 
G-)^{x,x') = —G^{x',xy. The quantum numbers of the 
quarkonium interpolating fields considered in this work 
are listed in Table |III| To optimize the overlap with the 
IS and IP ground states we use wave functions from 
a lattice potential model (see Appendix D of Ref. [26] ) 
in the construction of r(r) (as already mentioned in the 
previous section, the gauge configurations are fixed to 
Coulomb gauge). The heavy-quark mass used in the po- 
tential model is adjusted to match the mass used in the 
lattice QCD calculation. We use r(x) as the source for 
the quark propagator and a point source for the anti- 
quark. At the sink, the convolution in Eq. (13) is per- 



formed efficiently using fast Fourier transforms. 

In order to reach the high statistical accuracy needed 
to extract the small effects of the isospin charge den- 
sity, we compute quarkonium two-point functions for 64 
different spatial source locations (distributed on a cubic 
sub-lattice with spacing L/A) on each of the source time 
slices, and average over these source locations. Examples 
of free quarkonium two-point functions with the smear- 
ing technique discussed above are given in Figs. [3] and |4] 
Note that ground-state plateaus are reached already at 
a distance of ~ 0.4 fm in Euclidean time, which demon- 
strates the efficiency of the smearing technique used here. 



Name 


7^^^ 


r(r) 


rib 




0is(r) 


T 




0is(r) a J 


ht 


r+- 


0ip(r,i) 


XbO 




Ei<^ip(r,j)cri 


Xbi 


r++ 


Efe,; <^jki <Pip{r,k) ai 


Xb2 


rp + + 

^2 


(j>ip (r, j) (7k + (t>ip{r, k) Oj (with j / k) 



TABLE III: Smearing functions r(r) used in the quarkonium 
interpolating fields for the given representation of the cubic 
group, Ti, and values of parity, P, and charge-conjugation, C. 
The functions (j)is{Tc) and (j)ip{r,j) are eigenfunctions from a 
lattice potential model. 
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FIG. 5: Extracted inverse phase shifts for ?7i,-7r and T- 
TT scattering. Fitting the phase shift to pcot5(p)/m^ = 

h -T- , as shown by the shaded band, we can extract 

the scattering length shown by the point at p^/m^ = 0. 



D. Correlator ratios for energy shifts 

To investigate the effect of the medium on quarkonium 
propagation, we consider the correlators 

C{n;bb;t) - (O5,(0O„^+ (t)a|^(0)Ol^+ (0)) , (15) 

where (...) denotes path integration via the average 
over our ensembles of gauge configurations, and the in- 
terpolators C]j^+ and produce the quantum num- 
bers of n-pion and bb states as discussed in the pre- 
ceding sub-sections. States with the combined quan- 
tum numbers of the given quarkonium state {bb — 
r]b, T, hb, XbO, Xbi, Xm) and the n-pion system prop- 
agate in this correlator and naturally, the spectrum of 
this system is different from the sum of the spectra of n 
pions and of quarkonium because of interactions. At Eu- 
clidean times where only the ground state of the system 
is resolved (after excited states have decayed and before 
thermal states are manifest), this correlator will decay 
exponentially as 



C{n-bb;t)^Z^.-bbeM~KMt)^ 



(16) 
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where fi'jj.j^ is the ground-state energy of the combined 
system. 

To access the change in the quarkonium energy as a 
function of isospin density or chemical potential, we fur- 
ther construct the ratios 



R{n,bb;t) 



(17) 



Since the two terms in the denominator decay exponen- 
tially at large times as exp(— iJ^j^t) and exp(— i?„^+i) re- 
spectively, the ratio will behave as 



R{n;bb;t) 



(18) 



where AE^.-^^ = E^.-,. - E, 



bb 



— Et, is the quantity of 



central interest in our investigation. 



but they are indeed found to be small. The measured val- 



ues of the 5- wave phase shifts are given in Tables IV and 
\V\ while for the P-wave states we are unable to extract 
statistically meaningful results. Since the measured scat- 
tering momenta are small, it is possible to perform a fit 
to the effective-range expansion 



p cot (5(p) /m^r 



1 



(22) 



to extract the scattering length and effective range for 
these interactions. This extrapolation is shown in Fig. [5] 
and results in 111^^0.11^,7^ — 0.011(22) and mTr?",,^ ,r = 
—4(17) for the rji, state, and m^a-^.Tr = 0.017(37) and 
™7rT'T,7r = 2(29) in the case of the T. Clearly the inter- 
actions are extremely weak but are at the limits of our 
ability to resolve with current statistics. 



IV. QUARKONIUM-PION SCATTERING 

The quarkonium state in the presence of a single pion 
allows us to study the scattering phase shift of this two- 
body system using the finite- volume formalism developed 
by Liischer (28, 29J. The 5- wave quarkonium states we 
consider have angular momentum J = 0, 1 and define the 
total angular momentum of the entire system since the 
pion is spin-zero. Since the pion and bb states have differ- 
ent masses, the appropriate generalisation of the Liischer 
relation to asymmetric systems [3D] is required. We can 
define a scattering momentum p = |p| through the rela- 
tion 



asVWe + alMl + ^\aMVe + a'tM^ (19) 



atAE-r, + atM-r, + atM^ 



where Mj;^ = M^f'^^ is the kinetic mass of the bb state. 
The scattering momentum then determines the eigen- 
value equation 



P cot S-^.^{p) 



ttL \ 47r2 



S(x) = lim 



|n|<A 

E 



- An A 



n- 



(20) 
(21) 



that is satisfied by the bb-n scattering phase shift, 
^bb IT (p) ' scattering momentum. 

Since we have three different lattice volumes, we can 
extract the phase shift at multiple momenta. In Figure 
[5j we show the phase shifts that we extract for the J^fc-Tr 
and T-TT scattering channels. These interactions neces- 
sarily vanish in the chiral limit as the quarkonium states 
are chiral singlet objects. We therefore expect only small 
scattering phase shifts at the quark masses considered 
in our study. At the level of statistical precision we have 
achieved, we can barely resolve the phase shifts from zero. 



TABLE IV: The r/b-n phase shifts extracted using the Liischer 
method as described in the main text. 
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TABLE V: The T-tt phase shifts extracted using the Liischer 
method as described in the main text. 
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16^ 
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-0.0015(18) 


0.0037(56) 


0.0073(110) 


20^ 
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-0.0013(6) 


0.0063(34) 


0.0124(66) 


24^ 
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-0.0015(6) 


0.0128(57) 
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V. ISOSPIN DENSITY DEPENDENCE OF 
QUARKONIUM 

For larger isospin charge, we interpret the system of 
pions in terms of a medium of varying isospin charge den- 
sity once the ground state is reached. In the correlators 
C{n;bb;t), the quarkonium state exists in this medium, 
interacting with it. We consider first the 5- wave quarko- 
nium states as they are statistically better resolved than 
P-wave states. 



A. S'-wave states 

The correlators C(n, bb, t) are shown in Fig. [g] for 
bb = T at representative values of the isospin charge 
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and for a^m = 2.6 on the 20^ x 256 and 16^ x 128 en- 
sembles. The in-medium correlators on the 20^ x 256 en- 
semble exhibit a long region of Euclidean time in which 
they decay as a single exponential. This region overlaps 
with the regions in which the multi-pion correlators and 
the individual quarkonium correlators are saturated by 
their respective ground states. This gives us confidence 
that by considering the correlator ratios of Eq. (17 1 we 



can legitimately extract the quarkonium energy shifts in 
medium. On the ensembles with T = 128, thermal con- 
tamination is more significant and restricts the range of 
useful time-slices, particularly for large isospin charge. 

The correlator ratios, R{n,bb;t), discussed above, are 
shown for both T and r]b at a heavy quark mass agm ^ 
2.6 on the 20'^ x 256 ensemble for a range of different 
isospin charges, n = 6, 12, and 18, in Figs. [7] and |8] along 
with fits to time dependence using Eq. ( [l8| . Fits are per- 
formed over a range of times where both the individual 
multi-pion correlation functions and quarkonium corre- 
lation functions exhibit ground-state saturation and are 
free from thermal (backward propagating) state contam- 
ination. This is ensured by choosing the central fit range 
[^min, tmax] such that a fit ovcr the range [tmin-5, imax+5] 
has an acceptable quality of fit. On the 20^ x 256 en- 
semble, we choose tmin — 20 and t^^^ — 80, beyond 
which thermal contributions are apparent. Because ther- 
mal contributions are more significant for the ensembles 
with T = 128, we choose imax = 40 for these cases. Sta- 
tistical uncertainties are estimated using the bootstrap 
procedure. To estimate the systematic uncertainties of 
the fits, we calculate the standard deviation between 
the three energies extracted from fits with the ranges 

[^min 5; ^max 5]? [^min ? ^max] 5 and [tmin ~t~ 5,tniax ~t~ 5] 

on each bootstrap sample. The systematic uncertainty is 
then obtained as the average of this standard deviation 
over the bootstrap samples. 

The extracted energy shifts and uncertainties are 



shown in Table VI For larger values of n, the energy 
shifts become noisier and we limit our analysis to the 
range of isospin densities where a successful fit could be 
performed for a given ensemble. 

As a check of the methods of our study, we construct 
ratios in which we artificially remove the correlations be- 
tween the bb system and the many-pion state by evalu- 
ating Cif^{c)Cmv[c + 5c), where Cx{c) represents the 
correlation function for the quantity X measured on con- 
figuration c, and 5c is either a constant displacement or 
a random shift. In both cases, the removal of the corre- 
lation eliminates the signal for an energy shift. This is 
shown for the 77^ with n = 5 in Fig. [9] for random shifts, 
and the same qualitative effect is seen for all choices of 
the density and quarkonium state that are considered. 

To summarise the analysis of the correlator ratios for 
the S-wave quarkonium states, Fig. [TO] shows the isospin 
density dependence of the energy shift of the bb states 
from the free bb energy for both the T and 77b channels. 
Results are presented for the ranges of isospin charge 
density where a statistically meaningful extraction of the 



energy shift can be made. As can be seen, there is a 
significant negative energy shift for much of the range of 
isospin density that we have investigated. The magni- 
tude of this shift first increases as the isospin density is 
increased, before saturating at a value of about 1 MeV 
and finally decreasing to zero for large pi, albeit with 
increasing uncertainty. It is interesting to note that the 
saturation occurs at the point at which a marked change 
in the energy density of the many-pion system was ob- 
served in Ref. [6], and is likely caused by the changing 
nature of the screening medium at this point. The be- 
haviour at low densities is in line with the expectations of 
the potential model discussed earlier, but the saturation 
effect is not predicted (the model was developed using 
the measured shifts in the potential in the low density 
region, so this is not surprising). 

We have performed these calculations for all three en- 
sembles of configurations but have only been able to ac- 
cess a limited range of densities with the current sta- 
tistical precision. The results from all of the ensembles 
are consistent in the region in which they overlap. The 
16^ X 128 ensemble provides the largest density range. 

We also consider the shifts in the splitting between 
the ?7b and T energies in medium as a function of the 
density. We extract these shifts by calculating the cor- 
related differences between the individual energies using 
the bootstrap method. A summary of the isospin charge 
dependence of this splitting is shown in Figjll] There 
are hints that the T energy is shifted slightly more than 
the r/b energy by the presence of the medium. However, 
the smallest- volume data do not show this difference and 
it may not be significant. 



B. P-wave states 

We also analyse the lowest-energy P-wave quarkonium 
states, hh, XWi Xbi and Xb2, in medium. We find that we 
cannot resolve differences between the medium effects for 
these different states and so consider a spin average of 
their energies. It order to extract this average energy, we 
construct the product correlator 



^Xbo^Xbi^Xb 



(23) 



that decays with the average energy En^p = -^[SEn.ht + 

expression C^^ is 

short-hand notation for [C{n,bb;t)y ). We also consider 
the analogous S-wave spin average combination 



C's — C-^^Cy 



(24) 



These can be used to form correlator ratios, R{n, P; t) 
and R{n,S;t), analogous to Eq. (17). Since the P-wave 
quarkonium correlators are themselves statistically nois- 
ier than the S'-wave correlators (see Figs. [3]and|4]), the 
precision with which we can extract the P-wave energy 
shifts is reduced. In our fits to these correlators, we 
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FIG. 6: The correlators for the T in a medium corresponding to isospin charge n for n — 6, 12, and 18 are shown. Data are 
presented for asm = 2.6 on the 20'^ x 256 (upper) and 16^ x 128 (lower) ensembles. Correlators for the rjt in medium behave 
similarly. 



choose t„ 



15 and t„ 



35. Fig. [12] shows rep- where the potential shift is more significant. 



resentative correlators and correlator ratios for the spin 



average P— wave state, and Fig. 13 summarises the mea- 
sured shift, also showing the spin averaged S'-wave shift 
for comparison. The potential model expectation is that 
the P-wave shift will be larger than the S'-wave shift; to 
investigate this, we form correlated differences between 
the individual S- and P-wave shifts but the statistical 
uncertainty remains large and we are unable to resolve a 
difference. 



C. Heavy-quark mass dependence 



As discussed in Section III C we have performed cal- 
culations for four different values of the heavy-quark 
mass, a^m, ranging from the bottom-quark mass down 
to ~ 1.5 times the charm-quark mass. The analysis of 
the in-medium correlators and ratios is very similar for 
all masses and we do not present it in detail. To in- 
vestigate the variation of the energy shifts as a func- 
tion of the heavy-quark mass we compute AE"^ ^f^{asm) — 
AE^j^^^{as'm = 2.6) using the bootstrap method. Be- 
cause of correlations between the measurements for dif- 
ferent values of the heavy-quark mass, this provides a 
more statistically precise determination of the difference 
than would be evident from a naive comparison. Fig- 
ure [M] shows these energy differences for the different 
values of agm. It is apparent that the strength of the 
energy shift in both rji, and T increases as the heavy- 
quark mass decreases, in line with expectations from the 
potential model discussed above. Since the quarkonium 
states for lower heavy-quark masses are physically larger, 
they probe regions of larger quark-anti-quark separation 



VI. DISCUSSION 

Heavy-quark bound states provide an important probe 
of the properties of a medium and have been used in this 
work to investigate systems of large isospin charge den- 
sity created by many-pion correlators. Specifically, we 
have used lattice QCD to investigate how the presence 
of this medium modifies the NRQCD energies of various 
quarkonium states. Our calculations make use of ensem- 
bles of lattices with three different physical volumes at 
a single lattice spacing and at a single light quark mass 
corresponding to ~ 390 MeV. We have found a mea- 
surable decrease in the energy of both the rjt, and T states 
and in the spin-averaged P-wave energy. This decrease 
grows as the isospin charge increases before plateauing 
at a charge density at which Ref. [B] previously observed 
strongly non-monotonic behaviour of the energy density 
of the medium. The saturation of the energy shift pro- 
vides further support to the conjecture that a transition 
from a pion gas to a Bose-Einstein condensate of pious 
occurs at this point. In the region of low isospin den- 
sity where the energy shift is increasing, the quarkonium 
energy shifts are found to be qualitatively in agreement 
with our expectations from a potential model augmented 
with the hadronic screening effect found in Ref. [T5] . 

We have investigated how the observed energy shifts 
depend on the mass of the heavy quark-anti-quark pair, 
finding an enhanced effect for lighter masses. Given 
the phenomenological interest in J/^* suppression in 
medium, it will be interesting to investigate the analo- 
gous behaviour in the charmonium sector using alterna- 
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FIG. 7: The correlator ratios for the T in a medium corresponding to isospin charges n = 6, 12, 18. The shaded bands show 
the statistical uncertainties of fits of the form given in Eq. hm. Data are shown for asm = 2.6 on the 20^ x 256 (upper) and 
16'' X 128 (lower) ensembles. 




FIG. 8: The correlator ratios for the rjb in a medium corresponding to isospin charges n = 6, 12, 18. The shaded bands show 
the statistical uncertainties of fits of the form given in Eq. | |l8[ ). Data are shown for asm — 2.6 on the 20^ x 256 (upper) and 
16'^ X 128 (lower) ensembles. 
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FIG. 9: The ratio R{5,rjb;t) computed with and without the 
correct correlation between the r^b and many-pion system on 
the 20'^ X 256 ensemble, as discussed in the main text. The 
shift (difference from unity) is only apparent when correla- 
tions are included. 



tive formulations of the heavy-quark action more appro- 
priate for the charm quark. However, this is beyond the 
scope of the current work. 

A similar study of NRQCD quarkonium correlators in 
QC2D (two-colour QCD) at non-zero quark chemical po- 
tential was recently presented by Hands et al. in Ref . [31] . 
In contrast to QCD with three colours, in QC2D, the ad- 
dition of a quark chemical potential does not result in 
a complex action and numerical calculations can be per- 
formed efficiently [32H34] . In Ref. [6] it was pointed out 
that the phase structure of QCD at nonzero /x/ has an 
intriguing similarity to that of QC2D at nonzero quark 
chemical potential. It is apparent that the similarities 
persist to the case of quarkonium energy shifts in medium 
as an at least qualitatively similar dependence on the 
charge density/chemical potential is observed in the two- 
colour QCD case. Recent work [35H37] has probed the 
connections between different gauge theories with non- 
zero (isospin) chemical potentials and, as the extent of 
this similarity is surprising, this warrants further inves- 
tigation. 

Finally, by looking at quarkonium-pion correlation 
functions on three different volumes, we have extracted 
77h-7r and T-tt scattering phase shifts at low momentum, 
albeit with significant uncertainties. The interactions arc 
found to be weak in both cases as expected from chiral 
dynamics. 
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n 
















16^ X 128 


20^ X 256 


24^ X 128 


16^ X 128 


20^ X 256 


24^ X 128 


1 


-0.32(15)(25) 


-0.24(09)(04) 


-0.28(10)(04) 


-0.30(17)(31) 


-0.26(11)(05) 


-0.31(12)(01) 


2 


-0.76(21)(37) 


-0.48(12)(06) 


-0.46(14)(05) 


-0.79(23)(45) 


-0.57(16)(07) 


-0.52(15)(03) 


3 


-1.09(29)(47) 


-0.70(16)(07) 


-0.63(15)(07) 


-1.15(32)(58) 


-0.84(20)(09) 


-0.71(18)(05) 


4 


-1.32(37)(57) 


-0.85(18)(08) 


-0.76(20) (07) 


-1.40(41)(71) 


-1.03(23)(10) 


-0.85(23) (06) 


5 


-1.47(46) (66) 


-1.01(24)(08) 


-0.86(24) (05) 


-1.55(51)(84) 


-1.20(30)(08) 


-0.98(29)(06) 


6 


-1.56(56)(74) 


-1.11(34)(07) 


-0.97(33)(05) 


-1.64(62)(94) 


-1.30(37)(07) 


-1.10(39)(06) 


7 


-1.58(66)(77) 


-1.14(32)(06) 


-0.99(38)(09) 


-1.7(0.7)(1.0) 


-1.30(38)(05) 


-1.12(45)(08) 


8 


-1.55(77)(78) 


-1.20(44)(06) 


-0.99(47)(13) 


-1.6(0.8)(1.0) 


-1.33(54)(05) 


-1.10(57)(13) 


9 


-1.49(88)(76) 


-1.24(52)(05) 


-1.03(58)(21) 


-1.5(1.0)(1.1) 


-1.34(59)(06) 


-1.14(71)(17) 


10 


-1.4(1.0)(0.7) 


-1.22(52)(05) 


-0.92(65)(29) 


-1.5(1.1)(1.1) 


-1.33(59)(07) 


-1.01(80)(26) 


11 


-1.3(1. 1)(0.7) 


-1.05(66)(07) 


-0.84(72) (31) 


-1.4(1.3)(1.1) 


-1.11(75)(10) 


-0.91(87)(30) 


12 


-1.1(1.3)(0.8) 


-1.02(59)(07) 


-0.70(78) (43) 


-1.2(1.4)(1.1) 


-1.08(69)(11) 


-0.73(96)(43) 


13 


-1.0(1.5)(0.8) 


-1.03(68)(09) 


-0.5(1.0)(0.5) 


-1.1(1.7)(1.2) 


-1.10(78)(13) 


-0.5(1. 2)(0.5) 


14 


-0.8(1.7)(0.9) 


-1.01(74)(11) 


-0.4(1. 1)(0.6) 


-0.9(1.9)(1.3) 


-1.07(88)(17) 


-0.4(1.3)(0.6) 


15 


-0.6(1.9)(1.0) 


-0.86(83)(10) 


-0.2(1.3)(0.7) 


-0.7(2.1)(1.5) 


-0.90(99)(16) 


-0.1(1.5)(0.8) 


16 


-0.4(2. 1)(1. 2) 


-0.79(73)(14) 


0.0(1.4)(0.8) 


-0.6(2.3)(1.6) 


-0.81(90)(23) 


0.1(1.7)(0.9) 


17 


-0.3(2.3)(1.3) 


-0.87(87)(15) 


0.0(1.6)(0.8) 


-0.4(2.5)(1.7) 


-0.9(1. 1)(0.3) 


0.1(1.9)(0.9) 


18 


-0.1(2.5)(1.5) 


-0.83(73)(19) 


0.5(1.9)(1.0) 


-0.2(2.6)(1.9) 


-0.88(96)(36) 


0.6(2.3)(1.1) 


19 


0.1(2.6)(1.6) 


-0.88(77)(20) 


0.2(1.9)(0.9) 


0.0(2.8)(2.0) 


-1.0(1.0)(0.4) 


0.3(2.2)(1.0) 


20 


0.3(2.7)(1.7) 


-1.0(1. 1)(0.3) 


0.4(2.0)(1.0) 


0.2(2.9)(2.1) 


-1.2(1.4)(0.6) 


0.6(2.3)(1.1) 



TABLE VI: Fits to the S-wewe correlator ratios on the various ensembles for agin = 2.6. For each combination, we report: the 
mean and the statistical and systematic uncertainties. 
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